Pair Superfluidity of Three-Body Constrained Bosons in Two Dimensions 
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We examine the equilibrium properties of lattice bosons with attractive on-site interactions in the 
presence of a three-body hard-core constraint that stabilizes the system against collapse and gives 
rise to a dimer superfluid phase formed by virtual hopping processes of boson pairs. Employing 
quantum Monte Carlo simulations, the ground state phase diagram of this system on the square 
lattice is analyzed. In particular, we study the quantum phase transition between the atomic and 
dimer superfluid regime and analyze the nature of the superfluid-insulator transitions. Evidence 
is provided for the existence of a tricritical point along the saturation transition line, where the 
transition changes from being first-order to a continuous transition of the dilute bose gas of holes. 
The Berzinskii-Kosterlitz-Thouless transition from the dimer superfluid to the normal fluid is found 
to be consistent with an anomalous stiffness jump, as expected from the unbinding of half- vortices. 
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PACS numbers: 67.85.Hj,64.70.Tg,75.40.Mg 

Due to their remarkable versatility, cold-atom systems 
are ideally suited to realize quantum simulators for the 
physics of strong correlations |l| . A fascinating approach 
towards enhancing correlations in atomic systems steams 
from the fact that they can emerge via dissipative pro- 
cesses. In fact, dissipation-induced correlation effects 
were observed in an experiment with Feshbach molecules 
subject to strong inelastic collisions More recently, 
it was found that three-body loss processes for bosons in 
an optical lattice give rise to an effective Bose-Hubbard 
model description with a local three-body hardcore con- 
straint Q ■ Such a constraint stabilizes the system in the 
presence of strong attractive interactions, where dimer 
bound states proliferate. The effective lattice model de- 
scribing this situation is the Bose-Hubbard Hamiltonian 
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where t denotes the tunneling matrix element for near- 
est neighbor sites (ij) and U < an attractive on-site 
interaction Q. The filling n is controlled by varying 
the chemical potential /i, bi (b\) denote bosonic annihila- 
tion (creation) operators and i%i the number operator for 
bosons on lattice site i. In contrast to the usual Bose- 
Hubbard model, the Hilbcrt space is now restricted by 
the contraint (b\) 3 = to a maximum of two bosons 
on each lattice site. In another recent proposal Q, a 
similar effective lattice model of bosons with a three- 
body constraint was derived for spin-I atoms, which in 
addition includes an explicit correlated hopping term 

H' = —t' O^^j + h- c -) of the dimer bound states. 

The model in Eq. (1) exhibits an intriguing phase dia- 
gram, shown in Fig. 1 for the case of a two-dimensional 
square lattice. Recent analytical calculations I3L &J,\ 
and numerical works for the one-dimension case [3|,|5|, |6| 
exhibited that besides the trivially insulating phases at 
n = and n = 2, the system stabilizes two kinds of 
superfluid phases. The atomic superfluid (ASF) is char- 
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FIG. 1: (Color online) Left panel: Ground state phase dia- 
gram of the three-body constrained Bose-Hubbard model with 
attractive on-site interactions on the square lattice. Contin- 
uous (first-order) quantum phase transitions are shown by 
dashed (solid) lines. The cross indicates a tricritical point. 
The gray area contains the estimated DSF region. Right 
panel: Finite temperate phase diagram along n/\U\ = —0.5. 
The inset shows the temperature dependence of the super- 
fluid density p s at t/\U\ = 0.05 (DSF) for L = 10 and 20 
along with lines for normal 2T/n (dashed line) and anoma- 
lous 8T/tt (solid line) universal stiffness jumps. 



acterized by a finite atomic condensate with (pi) ^ 
and a finite superfluid response p s . For strong interac- 
tions \U\ ^> t however, a dimer (pair) superfluid phase 
(DSF) [8(] is stabilized, which is characterized by a van- 
ishing atomic condensate and (bi) = 0, but a finite dimer 
condensate density with an order parameter ({bi) 2 ) ^ 0. 
Such single-component DSF phases have been observed 
before in models with explicit correlated or pair hopping 
processes [1, |H,[l(J. In the DSF phase, the U(l) symme- 
try of the Hamiltonian is partially broken down to Z2, 
and at the DSF to ASF transition, this remaining Z2 
symmetry gets broken. With respect to this partial in- 
ternal symmetr y b reaking, DSF are thus related to spin 
nematic states [111 ]. Within Ginzburg-Landau theory, a 
Feshbach resonance term couples the ASF and DSF or- 
der parameter fields, which implies an effective U(l) x Z2 
symmetry similar as in the boson Feshbach resonance 



2 



problem [1 0, H El • From such analysis, the ASF-DSF 
transition was found to be Ising-like at unit filling, n = 1, 
and driven first-order by fluctuations via the Coleman- 
Weinberg mechanism 14 for n ^ 1 [j| Q . 

Here, we study the ground state and thermal phase dia- 
gram of the three-body constrained Bose-Hubbard model 
in Eq. (1) using quantum Monte Carlo (QMC) simu- 
lations. We focus on a two-dimensional square lattice 
geometry, on which true off-diagonal long-range order 
(ODLRO) within the ASF and DSF regimes emerges in 
the ground state. Besides establishing the presence of 
ODLRO, we assess the above mentioned theory for the 
ASF-DSF transition [H, 0] as well as a recent effective 
potential approach to the insulator- ASF quantum phase 
transitions [7| , which finds that they are driven first-order 
close to n = 1. We indeed obtain evidence for a tri- 
critical point along one of these transition lines. The 
thermal Berzinskii-Kosterlitz-Thouless (BKT) transition 
out of the DSF phase is found to be consistent with an 
anomalous jump in the su per fluid density, driven by the 
unbinding of half- vortices |l5| . 

Method- We employ a generalized directed loop al- 
gorithm in the stochastic series expansion (SSE) repre- 



sentation [16l |17I ] at finite temperatures T. The sim- 
ulations are performed on finite square lattices of lin- 
ear extend L (and N = L 2 sites), with periodic bound- 
ary conditions, such that the superfluid density is ob- 
tained from the winding number W fluctuations, p s = 
T(W 2 )/{2t) 0, [H. Since ODLRO in our system is 
forbidden at finite T [20| . we need to perform the sim- 
ulations at sufficiently low T to probe ground state 
properties of the finite system, as detailed below. Us- 
ing two kinds of directed loops, in which the worm 
heads carry either a single creation (annihilation) op- 
erator (b) or a pair operator (b^) 2 ((b) 2 ), provides 
direct access to the equal-time Green's function of the 
atoms Gi(i,j) = (b\bj) as well as the dimers G2(i,j) = 
(( b l) 2 (bj) 2 ) [SHil. The atomic and dimer condensate 
densities are obtained as C\ = ^/N 2 ^ i -Gi(i 1 j) and 
C*2 = l/N 2 J2ij G2{hj) respectively after extrapolations 
to the thermodynamic limit. A similar scheme with pair- 
worms was shown recently to be efficient for simulating 
two-component boson systems [22T |. Here, we find that 
accessing the dimer condensate density based on Gi still 
becomes problematic at the relevant low temperatures: 
Histograms of individual measurements for C2 (related to 
the lengths of the pair operator loops) exhibit fat-tailed 
distributions, i.e. the estimator of this quantity is dom- 
inated by rare events that make its sampling inefficient. 
A typical histogram within the DSF region is shown in 
Fig. 2. This behavior results form the fact that pairs of 
bosons proliferate at large \U\ near n = 1. A worm head 
carrying a bosonic pair operator performs an off-diagonal 
move (corresponding to the hopping of a dimer) only if it 
encounters a bond that shares a dimer (or an empty site, 
if the worm either carries annihilation or creation opcra- 
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FIG. 2: (Color online) Left panel: Finite size scaling of Ci and 
C2 for different values of t/\U\ from simulations at T/\U\ = 
0.01/L. Right panel: Extracting C2 from simulations of Hi 
and H2- Inset: Histogram P(C2) from simulations of H, with 
a power-law fit P(C2) oc (C2) a , a = —2.2 (dashed line) to 
the fat tail at t/\U\ = 0.045, p/\U\ = -0.5, L = 14 and 
T/\U\ = 0.002. 



tors) and a single atom. Such processes are thus strongly 
suppressed in the relevant parameter regime. Moreover, 
the fat tail of the histogram fits well to a Frechet dis- 
tribution. The exponent of the power-law decay in the 
tail is a > —3 (a w —2.2 for the histogram in Fig. 2), 
so that the variance does not exist, and the central limit 
theorem for the mean value does not hold. 

Hence, we resign to alternative estimates of the dimer 
condensate density in the DSF region. Wc employ two 
different means of adding correlated hopping terms to 
the original Hamiltonian that allow for efficient updates 
of the boson pairs. The results from both approaches 
agree within error bars, and with results based on G2 in 
cases, where the distribution of the loop lengths leads to 
a finite variance. In the first approach, we add to H the 
correlated hopping term H' , such that H\ = H+H J . The 
other approach couples C2 directly to the Hamiltonian, 

such that H 2 = H-h/NJ2i 3 (0j) 2 (^) 2 + h.c. + l) fea- 
tures correlated hopping terms between all sites of the 
system (the diagonal term allows the insertion of long- 
range vertices into the SSE operator string). The cou- 
pling h/N ensures an extensive energy. Besides the di- 
agonal update of the short-range terms in H2, we in- 
sert/remove long-range vertices using heat-bath proba- 
bilities (23|. We can now measure C'2 based on the esti- 
mator C 2 = (Y.b(H b ) 2 )/(Nt) 2 - (H h )/(Nh), where H h 
denotes all correlated hopping terms in H2 and Hi, the 
atomic kinetic energy term on bond b. This method re- 
mains robust down to very low values of h/\U\ ~ 10 -4 , so 
that we extract the condensate density of the model H by 
fitting to a low-degree polynomial (the data is found to be 
essentially linear up to t'/\U\, h/\U\ ~ 10 -2 ), perform- 
ing a bootstrap analysis for the error estimation. Such a 
scaling is shown in the right panel of Fig. 2, and extrap- 
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FIG. 3: (Color online) QMC data for Ci obtained for different 
system sizes along u/\U\ = -0.5 for T/\U\ = 0.01/L. The 
critical coupling (t/\U\) c ~ 0.054 is denoted by a vertical bold 
line. Right inset: data collapse based on three-dimensional 
Ising critical exponents. Left inset: Histogram P(W) of the 
winding number W taken in the ASF (t/\U\ = 0.06) and DSF 
(t/\U\ = 0.045) region at T/\U\ = 0.001. 



olations to the thermodynamic limit in the left panel. 

ASF and DSF phases - We concentrate on the fixed 
line p/\U\ = —0.5 (cf. Fig. 1) and vary t/\U\ in or- 
der to reach both the ASF and the DSF phases. The 
finite-temperature phase diagram along this line (cf. Fig. 
1) exhibits the temperature scale at which the system 
undergoes a BKT transition, below which superfluidity 
emerges. The transition temperature Tbkt was deter- 
mined from the finite size scaling of ps with system sizes 
up to L — 30, following Rcf. [24[. As seen from Fig. 1, 
the slope of the BKT- line changes near t/\U\ ~ 0.054. 
For smaller t/\U\, the system enters the DSF regime and 
the superflow is driven by virtual pair hopping processes. 
Correspondingly, we observe in the DSF region a strong 
even-odd effect in the winding number [|| (i.e. only even 
values are measured) , which is absent in the ASF regime, 
cf. the left inset in Fig. 3. This allows to locate the T > 
phase boundary separating the ASF and DSF phases (di- 
amonds in Fig. 1). Moreover, due to the nematic nature 
of the DSF, the BKT transition to this paired phase is 
driven by unbinding of half-vortices instead of the usual 
integer vortices as for e.g. the BKT transition to the 
ASF phase This leads to an anomalous univer- 

sal 8Tbkt/tt stiffness jump at Tbkt instead of 2Tbkt/tt 
as for the ASF phase [lfj. This anomaly is consistent 
with the finite size behavior of ps shown in the inset of 
Fig. 1, and was accounted for in the determination of 



Tbkt based on the scheme in Ref . [24[ . To assess the ex- 
tend of the ASF phase within the ground state, we mea- 
sured C\ for different system sizes at T/\U\ = 0.01/L, 
which we found necessary in order to access the finite 
system's ground state properties. The finite size data is 
shown in the main panel of Fig. 3. From extrapolations 
of C\ to the thermodynamic limit such as shown in Fig. 2 
(left panel), we find that below t/\U\ < 0.054 the atomic 



condensate density vanishes in the thermodynamic limit. 
Since ps remains finite below this value, we still expect 
the emergence of ODLRO from a dimer condensate Cy, 
beyond the ASF regime. In fact, the data of C2 shown in 
Fig. 2 extrapolates to finite values both within the ASF 
and the DSF region. 

ASF - DSF quantum phase transition - Upon fixing 
p/\U\ = —0.5, we cross the ASF-DSF quantum phase 
transition line slightly away from unit filling n = 1 (the 
density for t/\U\ = 0.054 is about 1.04), so it is ex- 
pected that the transition at fixed p/\U\ = —0.5 is first- 
order 0,0. However, as pointed out in Refs. @, the cor- 
relation length in the vicinity of the Ising critical point 
at n = 1 is expected to be large, diverging as (1 — n)" 6 . 
Close to n = 1, this severely exceeds the finite sizes ac- 
cessible to our simulations. Based on this scenario, the 
accessible finite systems could be still controlled by the 
adjacent Ising critical point at n = 1. We indeed ob- 
serve an approximate scaling in the data of C\ close to 
the transition point, described by the standard finite- 
size scaling ansatz Ci(L, r) = L~ 2l3 ^g(L 1 /"(T — r c )/r c ), 
where r = t/|£/|, and r c denotes the transition point. 
In the inset of Fig. 3, we show a corresponding data col- 
lapse of the available finite-system data for fJ,/\U\ — —0.5 
from the main panel, using three-dimensional Ising crit- 
ical exponents (3 = 0.3264(4) and v = 0.6298(5) [1| 
and T/\U\ = 0.01/L (for a dynamical critical exponent 
z = 1), giving (t/\U\) c = 0.054(1). Unfortunately, we 
are not able to discern whether the quantum phase tran- 
sition indeed is first-order or not (histograms of various 
observables did not exhibit any two-peak structures on 
the accessible system sizes L < 20 at these low tempera- 
tures). Controlling p as to fix n = 1 for all finite systems 
in order to directly address the n — 1 case also turned 
out unfeasible. Given the above estimate of (t/\U\) c , the 
mean- field value (t/\U\) c w 0.044 Q is found to under- 
estimate the stability of the DSF phase by about 20%. 
The calculations of Ref. Q, while accounting for quan- 
tum fluctuations effects, still show a similar deviation. 

ASF-insulator quantum phase transitions.- Next, we 
address the quantum phase transitions between the ASF 
and the (trivially) insulating phases at n = and n = 
2. While SF to insulator transitions are often continu- 
ous, the y c an be driven first-order for attractive interac- 
tions [7l.Ti2l.l26l]. We find from our QMC simulations that 
the transition from the empty system (n = 0) to the ASF 
is indeed first-order over an extended parameter regime. 
This is evident from robust discontinuities in both the 
filling n and superfluid density p Sl such as those shown 
in the left panel of Fig. 4 at p/\U\ = -1.2. From Ref. 0], 
the transition is expected to turn continuous beyond a 
tricritical point near p/\U\ w —0.9. Performing simula- 
tions down to p/\U\ = —1.5, we did not obtain evidence 
for a continuous transition, indicating strong effects of 
the Feshbach resonance coupling for p/\U\ < —0.5. The 
scenario of Ref. 0, including a tricritical point, is found 
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phase transitions in an extended three-body constrained 
boson lattice model, consistent with our findings. 
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FIG. 4: (Color online) Left Panel: Density n and superfluid 
density p 3 near the ASF-insulator transitions at (J,/\U\ = —1.2 
(T/\U\ = 0.008). Right Panel: Hole density (2 - n) at [i = 
and T/|Z7| = 0.08/L for various system sizes fitted to (2— n) = 
a(r — Tc) ln(6/(r — r c )) for the dilute bose gas transition at 
t c = 0.125. The inset shows the ratio (2 — n)/p 3 . 



to be realized for the transition from the ASF to the full 
system (n = 2), which we find to be continuous beyond 
fi/\U\ « -0.3, cf. Fig. 1. The right panel of Fig. 4 shows 
the filling n and the superfluid density p s along the line 
fx = 0, where no discontinuities are observed. To assess if 
this transition corresponds to that of a dilute bose gas of 
holes 27J , we analyze the behavior close to the transition 
point in more detail. In that case and for two dimensions, 
the density (2 — n) of the holes should exhibit a linear in- 
crease (2 — n) = a(r — r c ) ln(b/(r — r c )) with a logarithmic 
correction |28| in the vicinity of the transition point at 
t c = 0.125, while p s = (2 - n) in the dilute limit (HI- As 
seen from Fig. 4 (right panel), both observables fit well to 
these functional forms; the dynamical critical exponent 
equals z = 2 at this transition j27|, while our QMC tem- 
perature scales are sufficiently low to sample the finite 
system's ground state. 

Conclusion - The DSF in the considered boson lattice 
model is restricted to a narrow region of parameter space; 
the anomalous jump in the stiffness at the BKT transi- 
tion should however provide a direct experimental sig- 
nature of this phase. The ASF-insulator transitions are 
found in overall agreement with the effective potential ap- 
proach 0, which however underestimates the Feshbach 
resonance coupling effects at low filling. We established 
the presence of a tricritical point along the ASF to insu- 
lator transition line; determining its precise location and 
critical properties remain challenges for future studies. 
It will also be interesting to explore finite temperature 
transitions between the DSF and the ASF. 
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Note added.- After the completion of our simulations, 
we became aware of recent results 129| on the thermal 
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